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The random antiferromagnetic spin- 1/2 XX and XXZ chain is studied numerically for varying 
strength of the disorder, using exact diagonalization and stochastic series expansion methods. The 
spin-spin correlation function as well as the stiffness display a clear crossover from the pure behavior 
(no disorder) to the infinite randomness fixed point or random singlet behavior predicted by the the 
real space renormalization group. The crossover length scale is shown to diverge as £ ~ D -7 , where 
T> is the variance of the random bonds. Our estimates for the exponent 7 agrees well within the error 
bars with the one for the localization length exponent emerging within an analytical bosonization 
calculation. Exact diagonalization and stochastic series expansion results for the string correlation 
function are also presented. 

PACS numbers: 



I. INTRODUCTION 

Quantum spin chains exhibit a number of interesting 
features, especially at low temperature when quantum 
fluctuations are stronger than thermal ones. The anti- 
ferromagnetic (AF) Heisenberg model in one dimension 
(ID) has been extensively studied since the discovery in 
1931 of the Bethe Ansatai for the spin S = h chain. In 
ID, the AF XXZ model defined by the Hamiltonian 
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with J > 0, exhibits a gap-less excitation spectrum for 
A e [—1,1] for S = h (and more generally for half inte- 
ger spins 2 ), whereas a gap opens up in the spectrum for 
integer spins 3 ". In ID, the quantum fluctuations prevent 
the formation of long-range order— but the correlation 
length of the model [Eq.QJ] is infinite and a quasi-\ong- 
range order (QLRO) emerges, with power-law decaying 
spin-spin correlation functions in the ground state (GS) 
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where a = x, y or 2, 



is the GS expectation value, 



1 — j-i/ir, with 



and the critical exponent r\ x ^ 
[i = arccosA. 

If the AF exchange couplings are position-dependent, 
or more generally distributed randomly according to a 
probability distribution V(J), the situation changes dra- 
matically. Indeed, the spin-1 chain described by the 
random-exchange XXZ Hamiltonian 
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has lost the translation symmetry and rare events in the 
chain dominate the low energy physics^l Note that the 
energy scale is set to unity by choosing mean values of 
random couplings equal to 1. For instance, (i) the ran- 
dom planar exchange (RPE) model with J z (i) = l,Vi 
and J±(i) random; (ii) the random z-z exchange (RZE) 
model with J±(i) — l,Vi and J z (i) random; (hi) the 
random exchange (RE) XXZ antiferromagnet for which 
J±(i) = J z {i) and are all random numbers. 

For the AF XXZ spin-1 chain, it has been shown by 
Doty and Fished that disorder is relevant and that any 
amount of randomness destroys the QLRO and drives the 
system from a line of pure fixed points to an infinite ran- 
domness fixed point (IRFP)£. The situation for higher 
spins S > h is more complicated since it depends^ on the 
parity of 25 and some issues are still under debate >iS*ii. 
Regarding the thermodynamic properties of the random 
spin-1 XXX antiferromagnet, a real space renormaliza- 
tion group (RSRG) scheme, introduced first by Ma, Das- 
gupta and Hui£ lead to a number of analytical results. 
In particular, independent of the initial distribution "P( J) 
of couplings the low energy properties at the IRFP are 
characterized by a dynamical exponent z = 00 and a GS 
which consists of a tensorial product of randomly long- 
range coupled dimers, the so-called random- singlet phase 
(RSP) 7 . In such a phase, the disorder averaged spin-spin 
correlation function is dominated by strongly correlated 
pairs and is therefore slowly decreasing, as a power-law 



GS r VRSP ' 
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where rjnsp = 2 for all spin components a (• T '■ denotes 
the average over the disorder and the sites i). On the 
other hand, in the RSP the typical correlations decay 
faster (i.e. with a stretched exponential) than the aver- 
age correlations. These analytical predictions, that we 
will recall in greater detail in Section 2, have been tested 
numerically several times using different methods. First, 
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Lanczos exact diagonalization (ED) computations for the 
model (0 with RPE by S. Haas et ali£ for small chains 
(L m ax = 18) with the Hamiltonian J^J with RPE ap- 
peared to confirm the universal behavior predicted in 
Eq.@, at least for strong disorder. Using free-fermion 
ED, the RE XX model with the couplings uniformly dis- 
tribut ed ov er the interval [0, 1], has been investigated in 
Refs. I r ll ll for large systems (up to 1024 spins) and the 
critical behavior of the correlation functions was found 
to agree well with the analytical predictions^ More- 
over, a numerical RSRG calculation using again a uni- 
form coupling distribution has also confirmed theses pre- 
dictions for the RE XXX modetii and other features 
of the RSP have been checked at finite temperature by 
quantum monte carlo (QMC) simulations on finite chains 
(L max = 96) in RefE 

Whereas the universal behavior described by Eqs.QJ 
appears o be well established, the recent density matrix 
renormalization group (DMRG) calculations^ for chains 
(with free boundary conditions) defined by Eq.Q with 
RPE caused a debate^ which we intend to settle in this 
paper. Indeed, the conclusions of the DMRG simula- 
tions presented in Ref. on systems up to 400 spins, 
quite similar to a previous one using smaller systems, 21 
disagree with the IRFP scenario in so far as a dependence 
of the exponent rj upon A and the disorder strength was 
claimed. The RSRG scheme is expected to be asymptot- 
ically exact, but finite size (FS) effects cannot be negligi- 
ble, especially for weak disorder, i.e., far away from the 
IRFP. Indeed, one can show that the RG flow toward 
the IRFP is controlled by a crossover characterized by a 
length scale £ which is disorder dependent and diverges 
when the disorder strength is approaching zero. Such a 
behavior has already been mentioned in, 20 supported by 
exact diagonalization calculations for large system sizes 
(up to 4096 spins). 

In this paper, we intend to shed more light on this 
disorder induced crossover in random quantum spin-^ 
chains and demonstrate convincingly via numerical stud- 
ies of several related models defined by Eq.© the con- 
sistency of FS effects and the IRFP scenario. Such a 
crossover from pure to random critical behavior is very 
common in disordered systems and is always relevant, in 
experiments as well as in numerical studies, when the dis- 
order is not too strong and the length scales that can be 
explored are not too large. A good understanding of the 
order of magnitude of the crossover length scale, in par- 
ticular its scaling behavior in dependence of the disorder 
strength, is therefore necessary in order not to be mis- 
led by the mere appeara nce of the experimental and/or 
numerical data (c.f. Refs ll9l2"o|) . 

The paper is organized as follows. In Sec. II, we first 
recall the analytical predictions of the RSRG scheme. 
Then, using the bosonization study of the weakly disor- 
dered spin-i chain j£ we establish a disorder-dependent 
length scale which is the localization length of the re- 
lated problem of disordered particle in IDiS We will see 
that this length scale is closely related to the crossover 



1 I Z I O I t 

9 1 - 2 « 23 9 3, 9 (a) 

I 

1 2 3 4 

C ^ • (b) 

FIG. 1: Schematic picture of the decimation procedure for a 
4-spins cluster, (a) The strongest coupling, between spins 2 
and 3, is decimated — > (b) The spins 2 and 3 are frozen into 
a singlet and the spins 1 and 4 are coupled by an effective 
coupling Ji_4 (see the text). 



length scale that emerges in our calculations. In Sec. Ill, 
we present the free-fermions ED results for the spin-spin 
correlation functions for various disorder strengths for 
system sizes (with periodic boundary conditions) up to 
4096 sites. The crossover length scale, which emerges 
naturally from the data analysis, is studied as a func- 
tion of the strength of the disorder and compared with 
the localization length. The Ising part of the Hamilto- 
nian [Eq.@] has also been included and investigated via 
QMC simulations, using the stochastic series expansion 
(SSE) method. Section 4 is devoted to SSE calculations 
performed at very low temperatures for the RE XXX 
model and the RPE model. First a brief explanation of 
the method is given and some technical issues about equi- 
libration and GS convergence are discussed; then results 
for spin-spin and string correlation functions are shown. 
Finally in section 5, we give a summary and some con- 
cluding remarks. 

II. ANALYTICAL PREDICTIONS 

A. Real space renormalization group results 

The RSRG method, introduced originally by Ma, Das- 
gupta and Hu for the RE XXX spin-i chain 12 has been 
developed and studied in an exhaustive way by D. S. 
Fisher 7 for more general RE XXZ Hamiltonians. The ba- 
sic ingredient of this decimation procedure is a successive 
decrease of the energy scale via a successive decimation 
of the strongest couplings in the chain. Starting with a 
4-spin cluster (see Fig2J in which the exchange J2-3 is 
much larger than its neighbors J1-2 and J3-4, the spins 
2 and 3 are frozen into a singlet and the spins 1 and 4 
are coupled via an effective coupling 

t ^1-2^3-4 ,_s 

■A-4 = ^77 , (5) 

2J2-3 

calculated in second order perturbation expansion. 

In the limit of infinite system size, Fisher has demon- 
strated the existence of a fixed point for the the distribu- 
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tion of the effective couplings, independent of the initial 
distribution, which is given by 



Vo(J) oc J 



-1+8- 
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Such a distribution is very broad and guarantees that 
the RSRG procedure is asymptotically exact since the 
decimation produces weaker and weaker couplings. The 
IRFP, characterized by the distribution © , is attractive 
for any amount of randomness in the case of spins \ and 
the RSP, discussed in Section 1, describes the GS. At the 
critical point, the energy and length scales are related via 



lnA e ~ -vT 



(7) 



and as a consequence, the dynamical exponent z is in- 
finite. Concerning the correlation functions, as already 
mentioned in Section 1, the average and typical values be- 
have quite differently since rare events control the physics 
(see Eq. The average correlation function is domi- 
nated by long-range paired singlet and takes the following 
expression, independently of the direction (transverse or 
longitudinal) 

C avg (r) ex L_L. ( 8 ) 

Another quantity, which measure an hidden order is 
the string correlation function, defined in the GS at dis- 
tance r by 



S(r) 
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At the IRFP, the disorder averaged expectation value 
S avg (r) is expected to decrease as a power law, with a 
well-defined exponent^ 



Savg(r) = S(r) oc 



y.2—(f) 



(10) 



(p = being the golden mean. The RSRG method 

yields several other properties at the IRFP ; in particular 
finite temperature properties are well understood, but it 
is not the aim of our work. 

To give a brief conclusion for this part, we could say 
that this real space decimation procedure is analytically 
solvable for the random spin-i chain and is asymptoti- 
cally exact. Nevertheless, this theory does not quantify 
the FS effects, which are always present in numerical sim- 
ulations. Therefore, in order to give good interpretations 
of numerical results, the understanding of FS effects is 
crucial. This is what we are striving for here, using the 
bosonization treatment of the random chain. 



B. Bosonization of the random chain: Emergence 
of a disorder-dependent length scale 



In this part, we summarize previous results obtained 
using bosonization techniques AS The XXZ spin-i chain 



can be mapped using the Jordan- Wigner transformation 
(see Section 3) into a spinless interacting fcrmions prob- 
lem in ID. The low-energy excitations around the Fermi 
points can be considered in terms of bosonic fields and 
the resulting Hamiltonian describes a Luttinger liquid^ 
It is characterized by a set of A-dependent Luttinger liq- 
uid parameters which are the velocity of excitations u 
and the parameter K , given by 



u(A) = Z^, 
2 /j, 

K(A) 



2(tt - fx) 



(11) 



Several types of quenched randomness added to the pure 
XXZ model has been studied by Doty and Fisher in- 
where they found, for random perturbations that pre- 
serves the XY symmetry, a critical behavior which be- 
longs to the universality class of the Giamarchi-Schulz 
transition for ID bosons in a random potential^ Let us 
define the disorder strength V by 



V, 



(12) 



More precisely, for the RPE model 2?rpe = T > ±, for the 
RZE model, X> RZE = V z and for the RE XXZ model, 
since the randomness is isotropic £>re-xxz = = Dz- 
For a weak random perturbation added to the planar 
exchange, the renormalization under a change of length 
scale I = In L is 



dV 



= (3 - 2K)V. 



(13) 



Therefore, if K < 3/2 (i.e. — | < A < 1) the disorder 
is a relevant perturbation and the phase is the RSP. The 
renormalization flow toward the IRFP is controlled by a 
length scale which emerges from Ea. (|13|l : 



f*(X>) ~ v- 



(14) 



For the RE XXX model, the random perturbation added 
to the operator S*S* +1 is marginally irrelevant and there- 
fore the exponent 3 _l )K — \ for A = 1 is expected to 
have small logarithmic corrections^. 

The length scale £* is called the localization length 
since in the fermionic language, the transition at T> > 
is a localization transitioni22iSi Such a metal-insulator 
transition driven by the disorder is characterized for in- 
stance, by the vanishing of the zero temperature Drude 
weight (also called the charge stiffness or the spin stiffness 
in the case of the spin- 1/2 XXZ model) V2? > in the 
thermodynamic limit. Previous numerical studies have 
checked this effect on transport properties in the case of 
interacting spinless fermions in a random potential us- 
ing ED2£ or equivalently in the case of the XXZ chain in 
a random magnetic fieldi22*22i For the simpler model of 
non-interacting fermions with random hopping, mapped 
into the RE XX spin chain, very large scale numerical 
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simulations have been carried out on systems up to 2048 
sites' Using the scaling law for the spin stiffness 



Ps (L,V)=g(L/C(V)), 



(15) 



the localization length has been precisely studied and 
agrees perfectly with Eci. ((T4fl for weak disorder (see Sec- 
tion 3 C and Figs l5l6(l . Our purpose here is to study 
crossover effects for various ID spin-^ disordered mod- 
els. Regarding the low-energy effective theory predicted 
by some bosonization calculations, there is only one rel- 
evant length scale which emerges from it, i.e., the lo- 
calization length £* (T>) . Based on numerical calculations 
performed over FS clusters for various disorder strengths, 
the next sections are dedicated to the study of the dis- 
order dependence of the crossover length scale and its 
comparison with the localization length. 



III. EXACT DIAGONALIZATION STUDY AT 
THE XX POINT 

A. Free fermions representation 
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FIG. 2: Averaged correlation function Ca Vg (L/2) as a func- 
tion of the system size L on a log-log scale for W = 
0, 0.25, 0.5, 0.625, 0.75, 1.0 (from top to bottom). The data are 
averaged over 50000 samples for L < 1024, 3000 for L = 2048 
and 500 for L — 4096, the statistical errors are smaller than 
the symbol sizes. The data for the pure system (W = 0) follow 
C x (L/2) tx L~ 1/2 , the full line with slope -2 is the expected 
asymptotic behavior according to the IRFP scenario. 



Let us consider the ID XX spin \ model with ran- 
dom exchange couplings J±(i). This quantum problem 
is governed by the following lattice Hamiltonian 



i=l 



(16) 



We impose periodic boundary conditions; Sl+i = Si. 
It's well known that this spin problem can be mapped 
into a free spinless fermions model via the Jordan- 
Wigner transformation : Sj = 1/2 — rij, and 5+ = 

Ej-i 
i=1 c sa tisfy fermionic commutation re- 

lations, {c\,Cj} — Sij, and rij = c^Cj is the number of 
fermions (spin down) at the j-site. The Hamiltonian can 
then be written as 
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(17) 



where h.c is the hermitian conjugate and AT = Y]f—i is 
the number of fermions in the system. In the non random 
case, the solution of the problem via a Fourier transfor- 
mation is trivia!^ because of the translational invariance. 
But in the random system, this symmetry is broken and 
we have to solve numerically a random matrix problem. 
The way to obtain the correlation functions is straightfor- 
ward and has already been explained in several previous 
works it amounts to a numerical calculation of 
the eigenvectors of a L x L band matrix and then the 
evaluation of a (r — 1) x (r — 1) (resp. 2x2) determinant 



in order to compute the transverse (resp. longitudinal) 
spin-spin correlation function at distance r C x (r) (resp. 
C z (r)). We can note that in the same way, the string 
correlation functions can also be obtained^ 



B. Numerical results for the spin-spin correlation 
functions: crossover effects 



In order to study the crossover as a function of the 
disorder strength, we have chosen the following In- 
dependent flat bond distribution 



P{J±) 



otherwise 



W, 1 + W] 



(18) 



The disorder strength, defined by Eq.JT^J, is V = W 2 /3 
and we define S as the variance of the random variable 
In Jj_(i) by 



<5 2 = (lnJ_L(z)) 2 
which is related to W according to 



(lnJx(i) 



1 - 



1 - W 2 
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(20) 
D. In 



We note that for weak disorder W <C 1, 5 ' 
order to reduce statistical errors and boundary effects we 
have used the PBC and computed the bulk correlation 
function in the transverse direction at mid-chain 



(21) 
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for several system sizes (L = 2 q , q = 1, 12) and disor- 
der strengths (W = 0.25, 0.5, 0.625, 0.75, 1). The data 
for C x (j) were computed for each individual sample us- 
ing standard routines and then averaged over the disor- 
der. The number of disorder configurations was more 
than 5 • 10 4 for L < 1024 and at least 500 for the largest 
size and weakest randomness (L — 4096 W = 0.25). 
In FigEl w e show the average bulk correlation function 
Cavg(%) — C x (^) for different disorder strengths. We 
observe that for small system sizes the slope of C x {\) 
versus L in a log-log plot is much smaller than 2, the 
value that one would expect form the IRFP scenario. But 
when L increases one observes a crossover, as reported in 
Ref. |20j, from an apparently non-universal behavior with 
a FF-dependent power law exponent r](W) for small sizes 
to a universal behavior with C^ vg (j) ~ L~ 2 for L — > oo, 
as predicted by the RSRGi£ Such a behavior suggests the 
existence of a disorder-dependent crossover length scale 
£ which controls the crossover from the pure (instable) 
fixed point to the IRFP which is attractive, even for weak 
disorder. Defining the dimcnsionless parameter x = L/£, 
one can identify three different regimes: 

(i) For x -C 1, the critical behavior of the pure sys- 
tem (J±(i) — constant) is dominant, with an exponent 
r,(W) = 1/2. 

(ii) For x 3> 1, we are in the asymptotic regime where 
the predictions of the RSRG are recovered, in particular 
V(W) = r) RSP = 2. 

(iii) For x ~ 1 we are in the crossover regime with a W- 
and L-dependent effective (FS) exponent r](W). 
Consequently, we expect the following scaling form: 
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FIG. 3: Scaling plot according to Eq. (1221 of the data shown 
in FigHwith £ = 600, 140, 88, 54, 20 for W = 0.25, 0.5, 0.625, 
0.75 and 1.0, respectively. The symbols are identical to Figl2l 
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where the scaling function c{x) is constant in the regime 
(i), and c(x) — > x~ 3 / 2 in the regime (ii). In Fig|3 the 
scaling plot following Ea. l22|) is shown for the data of 
Fig0 f (W = 1) has been chosen such that the crossover 
region (iii) is centered around x ~ 1 and the other esti- 
mates have been adjusted in order to give the best data 
collapse. 



FIG. 4: Disorder dependence of the crossover length scale £ 
of the random XX chain. The full squares are the numerical 
estimates from the data collapse in FigE] and full lines are 
fits, (a) In function of the disorder parameter T>, a power law 
with an exponent —1.048 fits the data only for weak disorder 
whereas in (b), a fit £(<5) ~ <5 -1,8 works for the entire range 
of disorder strength studied here. 



C. The crossover length scale as a localization 
length 

In this part, the dependence of the crossover length 
scale on the disorder strength is studied. A compari- 
son with the localization length £* , calculated using the 
spin stiffness of the RE XX chain, is also presented. 
Fig El shows a plot of £ vs the disorder parameters V 
and 5. As expected one can observe a singular behavior 
for T> or 5 — > 0. More precisely, we observe in FigQfa) 
that for sufficiently weak disorder (typically for I? < 0.1), 
the crossover length scale is well fitted by a power-law : 
£(£>) oc T>~ 7 with an exponent 7 = 1±0.1, in good agree- 
ment with the localization length exponent predicted by 
Ea.lfTlfl which gives 3 J L 2K = 1 at the XX point. For 



stronger disorder, a deviation from the power-law is ob- 
served. On the other hand, £(<5) shown in Fig^b), can 
be fitted by a power-law £(S) cx 5"*, with $ = 1.8 ± 0.2 
for the whole range of randomness studied here. 

It is instructive to compare the crossover length scale 
£ with the localization length £*, extracted from the 
numerical calculation of the spin stiffness of the RE 
XX chain (for more details about this calculation, see 
Ref. Isoh . While the transport properties of random spin 
chains arc not the purpose of this paper^ we mention 
here some results that two of us obtained by ED per- 
formed on the RE XX chain*2 The spin stiffness p s which 
measures the magnetization transport along the ring is 
calculated in the GS as the second derivative of the GS 



6 



1 Q. 



0.1 - 



T 



— p s =l/7t pure case 

* W=0.025 ^'=7875 
■ W=0.05 % =2135 

W=0.075 (' =882 

* W=0.1 % =500 
W=0.125 ^*=332 
W=0.15 ^'=228 
W=0.175 E,*=168 

+ W=0.2 I* =133 
W=0.225 |*=100 

* W=0.25 % =81 



W=0.3 
■ W=0.35 
♦ W=0.4 

W=0.45 
■* W=0.5 

(L/i;*)' ' 5 
L 



£,'=55 
c/=44 
i"=30.5 

k* =23 
|*=17.5 
IRFP 




0.01 



10000 



100 



FIG. 5: Inverse logarithm of the disorder averaged spin stiff- 
ness plotted for several box sizes W specified on the plot. 
All the curves are collapsed since a rescaling of the x-axis has 
been done, providing an universal curve as a function of L/£* . 
The W-dependent localization length £* has been calculated 
for each disorder strength, as indicated on the plot, in order 
to give the best data collapse. The full line stands for the 
pure case and the dotted one is for the IRFP behavior. 
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FIG. 6: Disorder dependence of the localization length 
of the random XX chain calculated using the scaling of the 
stiffness [Ea. l|15|l ]. (a) In function of the disorder parameter 
T>, the expected power- law Eg. 1141 with an exponent equal to 
— 1 is in perfect agreement with numerical data (open circles) 
which can be fitted for weak disorder, with an exponent equal 
to — 1±0.01. (b) In function of the other disorder parameter S, 
the numerical data (open circles) are perfectly described by a 
power-law, for the entire range of disorder, with an exponent 
equal to -2 ± 0.02. 



energy per site with respect to a twist angle tp applied 
at the boundaries using the so-called twisted boundary 
conditions^ and taking the limit p — > 0. For the same 
model [Eg . H16f l] studied in this section and for systems 
sizes going from 8 to 2048 sites, p s has been calculated 
by ED techniques for several disorder strengths (from 
W = 0.025 to W = 1) and averaged over a very large 
number of samples (from 10 5 for the smallest sizes to 500 
for the largest). 

The stiffness p has dimension of inverse (length d ~ 2 x 
£ T ), where £ r is the correlation length in the imaginary 
time direction^ In our case £ T ~ exp(A^ 1 ^ 2 ), which is 
one manifestation of the IRFP that dominates the critical 
behavior of the random XX chain, and £ = L for a finite 
system at criticality we expect p to scale as In ps(L) ~ 

— y/L. Combining this with Ea. ()15fl we show in FigEla 
scaling plot of — (In <7(L/£*)) _1 which displays the same 
features as FigEI Indeed, for L <C £*, the pure behavior 
is observed with a stiffness ps — and for L 3> 
£*, the IRFP behavior is recovered with In g (£/£*) ~ 

— (L/t;*) - 5 ; the regime where L ~ £ * being a crossover 
regime. 

The localization length £*(W) has been estimated for 
different values of the disorder strength (note that the 
computational demand for calculating the stiffness is sub- 
stantially smaller than the one for the correlation func- 
tion^, for which reason we could compute more data 
points) and is shown on the FigElversus the disorder pa- 
rameter. We see clearly that the behavior of the crossover 
length £ as a function of the disorder strength (see Fig0J) 



is exactly analogous to the one of the localization length 
£*. Indeed, for V <C 1, the bosonization result Eq.(nj 
agrees with numerical results, as shown in FigEJa), and 
for stronger disorder we observe the same deviation as in 
Fig- Etb) gives us the confirmation that for strong 
disorder the equation l|14l) has to be replaced by 



£*(($) ~ J ■ ( 23 ) 



Since for weak disorder 8 ~ Vv, we expect: <I? = 3 \ K 
which works perfectly for the entire range of disorder 
considered here, as shown in FigEIb). 

Let us summarize our results that we obtained so far 
for the RE XX chain. With ED calculations we stud- 
ied the crossover that controls the renormalization flow 
starting from a system with a finite disorder to an in- 
finite disorder fixed point. As predicted by RSRG and 
bosonization calculations, the IRFP is attractive for any 
amount of initial disorder and the crossover length scale 
£ is well described by a power-law, diverging like 2?~ 7 . 
Moreover, the exponent 7 has been identified to be iden- 
tical with the localization length exponent occurring in 
£* (D) ~ V~ , While the parameter V is suitable to 
quantify the divergence near 0, we have found the pa- 
rameter (5, Eq. (fHty\ to be a better candidate to describe 
localization and/or crossover behaviors for any strength 
of randomness : £(<5) — £*(<5) cx <5 - * with $ = 3 J 2K . 
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IV. QUANTUM MONTE CARLO STUDY 

A. The SSE method 

The SSE QMC method has been described as a loop al- 
gorithm in detail by one of us in^l More recently the con- 
cept of directed loop has been developedS&iSiiP. and the 
efficiency of such an algorithm has been demonstrated for 
several models, in particular for the XXZ model, defined 
by Eq.Q. We start from the general random-exchange 
XXZ Hamiltonian JSJ that we can rewrite as a sum over 
diagonal and off-diagonal operators 



n 



XXZ 
random 



J2[Mb)H hb - J±(b)H 2 . b 



6=1 



where b denotes a bound connecting a pair of interacting 
spins (i(b),j(b)). 

H hb = C - AS? (b) S? (b) (24) 
is the diagonal part and the off-diagonal part is given by 



Hi,. 



i S ?(b) S j( b ) 



i(b) 



(25) 



in the basis {\a >} = {\Sf, 5|, >}. The constant 
C which has been added to the diagonal part ensures 
that all non-vanishing matrix elements are positive. The 
SSE algorithm is based on Taylor expanding the partition 



function Z = Tr{e 



-put 



■ } up to a cutoff M. which is 



adapted during the simulations in order to ensure that 
all the elements of order higher than M. in the expansion 
do not contribute. So 



f3 n (M-n)\ 
Ml 



M 



alJlJMHa,^ a), (26) 



i=i 



Metropolis acceptance probabilities 



Pi 



[0,0] p ^[l,6] p 



^l,6] p — [0,0], 



J z (b)Lp(a(p) Hi, b a(p) 



M — n 



M 



1 



J z (b)L0(a(p) H ltb a(p) 



,(28) 



(29) 



During the "propagation" fx 
"propagated" state 



'om 



p = 1 to p = M, the 



(30) 



is used and the number of non-unit operators n can varies 
at each index p. The following step is the off-diagonal up- 
date, also called the loop update, carried out at n fixed. 
Its purpose is to substitute [1, bi) p <-> [2, bi] p in a cluster- 
type update, i.e., with the operators forming closed loops. 
Such a construction has already been discussed in detail 
elsewhere^. A very efficient directed loop implementa- 
tion can be used and for A € [0, 1] it has been shown 
that during the construction of the loop, back-tracking 
processes can be avoided. At the SU(2) AF point, the 
algorithm is deterministic because we can build all the 
loops in a unique way. So, for A = 1, all the loops 
are updated independently of each other with probabil- 
ity 1/2. For A 7^ 1 the construction of the loop depends 
on some well defined probabilities^ at each time a non 
unit operator is encountered in the loop building. 
One MC step is consists of diagonal updates at all possi- 
ble locations in the index sequence, followed by a number 
of loop updates (the number adjusted so that the average 
number of operators changes is comparable to the total 
number of operators). Before starting the measurement 
of physical observables, one has to perform equilibration 
steps, notably necessary to adapt the cutoff A4. 



where Sm denotes a sequence of operator indices 



B. Convergences issues 



Sm = [a>i,h], [02,62], ...[om^m] (27) 

with a, = 1,2 corresponds to the type of operator (diag- 
onal or not) and 6j = 1,2, ...L is the bond index. Note 
that ,h(b) = J z (b) and J 2 (b) = J±(b). A Monte Carlo 
configuration is therefore defined by a state \a) and a 
sequence Sm- Of course, a given operator string does 
not contain M. operators of type 1 or 2, but only n; so 
in order to keep constant the size of Smj M. — n unit 
operators -ffo.o = 1 have been inserted in the string, tak- 
ing into account all the possible ways of insertions. The 
starting point of a simulation is given by a random ini- 
tial state I a) and an operator string containing M. unit 
operators [0, 0]i, [0, 0]m- The first step is the diagonal 
update which consists in exchanging unit and diagonal 
operators at each position p [0, 0] p «-> [1, bf\p in Sm with 



The precise determination of physical observables us- 
ing QMC suffers obviously from statistical errors since 
the number of MC steps is finite. As we deal with disor- 
dered spin chains, the sample to sample variation is an- 
other source of errors. Moreover, the calculation of GS 
expectation values for a system close to an IRFP, where 
FS gap scale like lnA e ~ >/L, requires a very carefully 
numerical treatment. In order to avoid finite tempera- 
ture effects and to ensure that we measure observable 
in the GS, we use the /3-doubling scheme, developed in 
Ref. ^3 and then used in Refs. lllll42L Such a scheme is a 
very powerful tool because it allows to reach extremely 
low temperatures rather rapidly and reduces considerably 
equilibration times in the MC simulation. The proce- 
dure is quite simple to implement and its basic ingredient 
consists in carrying out simulations at successive inverse 
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FIG. 7: Test for the convergence of the disorder averaged 
longitudinal correlation function calculated for a the RE XXX 
chain at W — 0.5 with 16 sites. Results, averaged over 10 3 
samples, are for different number of MC steps (N e q,Nm) as 
shown on the plot. The /3-doubling scheme has been used 
with inverse temperatures /3 n = 2", with n used here for the 
x-axis. 




P 

FIG. 8: Test for the GS convergence of C7J„ g (~), defined by 
Eq.J3U, vs the inverse temperature /3. SSE calculations per- 
formed on the RE XXX model for W — 0.6 using the /?- 
doubling scheme with (N eq ,N m ) = (50, 100). Averaging has 
been done over 10 3 different samples and the results are shown 
for the 4 larger sizes L — 32, 48, 64, and 96. In the inset, the 
GS inverse temperature /3gs (see the text for its definition) is 
plotted in a log scale vs the the square root of system sizes. 
A linear fit is represented by the full line. 



temperatures j3 n = 2™, n — 0,1, ...,n mQX . Starting with 
a given sample at n = we perform a small number of 
equilibration steps N eq followed by N m = 2N eq measure- 
ment steps. At the end of the measurement process, (3 
is doubled (i.e. n — ► n + 1) and in order to start with 
an "almost equilibrated" MC configuration, the starting 
sequence used is the previous Sj^i doubled, i.e., 

S2M = [ai,h], ...[a M ,bM}[aM,bM], [ai,h]. (31) 

Such a scheme becomes very efficient at low temperature 
and for disordered systems, in which very small correla- 
tions may develop. It is for the moment the most efficient 
technique available to cancel remaining temperature ef- 
fects although a zero-temperature SSE algorithm might 
be developed soon4£ The next point concerns the num- 
ber of equilibration and measurement steps that we have 
to perform. It is illustrated for an L = 16 RE XXX chain 
with random bonds distributed according to Eq. (|18fl with 
W = 0.5 in Fig[7| Here the disorder averaged mid-chain 
longitudinal correlation function 



L 2 v— 

Cl V g{-^) = £ ^ < SiS* + L >GS (32) 

i=l 

is plotted for different values of (N eq ,N m ). The averag- 
ing is done over 10 3 independent samples and we observe 
that when the temperature becomes low enough, even for 
a couple (N eq: N m ) quite small, averaged values do not 
depend on the number of MC steps. As already men- 
tioned injii^ii^ we conclude that the sample to sample 



variation produces larger error bars than statistical er- 
rors, even for a number of measurement steps < 100, and 
in the following we will use the /3-doubling scheme with 
(N eq , N m ) — (50, 100) and a sufficiently large number of 
samples (> 10 3 ). 



In order to make reliable predictions for the GS, very 
large f3 have to be reached. This is illustrated for the RE 
XXX model with disorder strength W = 0.6 in Fig|8J 
where C^ vg (^) is plotted vs (3 for different chain sizes L. 
We consider that the GS expectation value is obtained 
when there are no statistically significant differences be- 
tween the results for (3 max — 2™ max and j3 — 2™ m< " _1 . 
More precisely, our GS convergence criterion is the fol- 
lowing : the GS is considered reached if the expectation 
value is 98% of the saturation value. Note that using 
such a criterion, we can define a system size dependent 
temperature scale below which the thermal expectation 
values are indistinguishable from GS expectation values: 
Pas = 2"" lQ *~ 1 ± 2 nmax ~ 2 and as shown in the inset of 
FigEl we obtain for this quantity a FS scaling of the form 
In pas ~ for W — 0.6. Note that we have checked 
the validity of this scaling for all disorder strengths con- 
sidered here. Such a scaling is not surprising since the 
FS gap also obeys to a similar law Eq.Q. 
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FIG. 9: Averaged longitudinal correlation function C z (L/2) 
for the random XXX model as a function of the system size 
L on a log-log scale for W = 0,0.25,0.5,0.6,0.8,0.9,1.0 and 
5 = 2 (top to bottom). The data, computed in the GS using 
SSE method and /3-doubling scheme, are averaged over more 
than 1000 samples. The data for the pure system (W = 0, 
open circles) follow C z (L/2) oc \/m(L)/L, the dashed line 
with slope —2 is the expected asymptotic behavior according 
to the IRFP scenario. 
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FIG. 10: Scaling plot according to Eq. 1351 for the data of 
the FigEHwith £ = 87, 37, 28, 16.5, 13, 10 for W = 0.25, 0.5, 
0.6, 0.8, 0.9 and 1.0, respectively. The full line stands for the 
pure behavior and the dashed line is the expected asymptotic 
behavior according to the IRFP scenario. 



case the exponent in Eq.|J2J is r\ z = 1, but logarithmic 
corrections have to be taken into account^ 



C. Spin-spin correlation functions of the 
random-exchange XXX model 

After these careful checks of equilibration and tem- 
perature effects in our simulations, we can analyze the 
SSE results obtained for the disorder averaged longitu- 
dinal spin-spin correlation C z vg . In order to extract the 
bulk value, we compute this quantity at mid-chain and 
perform the averages along the chains and over random 
samples, according to Eq. • We consider in the fol- 
lowing the RE XXX Hamiltonian 

L 

1~(-RE X = [^(*)(^f Sf+l + Si^i+l + Si^i+l > (33) 
i=l 

with J(i) random AF couplings taken from the In- 
dependent distribution Ea. lfTHl) . We have also used the 
more singular distribution V{ J) = 5J~ 1+S if J < 1 and 
otherwise, with 5 = 2. Such a distribution is, a priori, 
closer to the IRFP and therefore we expect the asymp- 
totic behavior C z vg (^) ~ L~ 2 to become visible already 
for not too large system sizes. Indeed, this is what we 
can see in FigEl where C z vg (^) is plotted versus L for 
different disorder strengths. The crossover phenomena, 
already mentioned for the RE XX case, is also clearly 
visible but from 16 sites the RSP behavior is recovered 
for the 5 = 2 case. For weaker disorder, the asymp- 
totic behavior is visible only for larger distances and an 
analysis analogous to the one we have performed for the 
XX chain is necessary in order to extract a disorder- 
dependent crossover length scale £. In the pure XXX 



C 2 (r) cx (-l) r 



(34) 



with which our numerical data for W = agree (see 
Fig©. 

As in the XX case we expect a disorder dependent 
length scale £ to govern the crossover from pure XXX 
behavior for L <C £ to the asymptotic RSP behavior vis- 
ible for L < £. Then C z {L/2) should obey the following 
scaling form: 



C z (— 

avg \ n 



L 



(35) 



with c(x) a scaling function that is constant in the 
pure regime (x 1) and proportional to (x ln 1 ^ 2 a;)™ 1 
for x 1 in order to reproduce the IRFP behavior 
OavAj) oc L- 2 for L > f. In FigUHl the scaling 
Eq. |3o|) is shown for the data of FigOU The W-dependent 
crossover lengths scale £ was chosen for each value for W 
individually to obtain the best data collapse. ^(W = 1) 
has been chosen such that the crossover regime is cen- 
tered around x ~ 1 (i.e. when the system size is of the 
same order of magnitude as the crossover length scale 
£). In comparison with the XX results (see section 3), 
the asymptotic behavior C z vg ~ L~ 2 sets in already for 
smaller system sizes. This observation is compatible with 
the fact that that the disorder dependent length scale de- 
fined in Ea. l|14l) diverges much slower at the XXX point 
(txxx oc V- 1 / 2 ), than at the XX point (£ X x ocP" 1 ). 

The disorder dependence of the crossover length scale 
of the RE XXX model is shown in FigEl For D -> 
it diverges with a power law and for small disorder 
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FIG. 11: Disorder dependence of the crossover length scale 
£ of the random XXX chain. The full squares are the nu- 
merical estimates from the data collapse in Fig llOl (a) In 
function of the disorder parameter T>, the power law fit 
£ ~ 5^ / °' 61±0 ' 1 works only for weak disorder whereas in (b), 
the fit £ ~ 8^" 16±0 ' 2 works for the entire range of disorder 
strength studied here. 



strengths, we can fit the data well by £(X>) ~ £>-o.6±o.i 
[see fig HT^ a)]. As a function of 5 the fit £(5) ~ <J- x - 2±0 - 2 
is working in the whole range of disorder strengths [see 

figHHb)]. 

The agreement of our numerical estimate of the expo- 
nent governing the divergence of the crossover lengths 
with the bosonization prediction for the localization 
length (0.6 ± 0.1 versus 0.5) is not as good as in the 
XX case but still acceptable within the error bars. These 
minor deviations might be due to small logarithmic cor- 
rections to formula (|14|) . This is expected since the 
bosonization approach gives predictions for the RPE 
model, whereas the RE case considered here is only qual- 
itatively similar because the randomness added to the 
Ising term is marginally irrelevant^ 



D. String correlation function 

The string correlation, Eq. @ , was introduced to mea- 
sure hidden order in in integer spin chains where the or- 
dinary spin-spin correlations vanish exponentially. In the 
RS phase the decay of the string correlation is expected 
to be described by a power law [see Eq.lJTU))]. with a de- 
cay exponent of ij ~ 0.382. It has been shown befor o 11 ! 14 
that the string correlation converges particularly quickly 
to the expected behavior. 

In this section we begin by demonstrating yet another 
crossover effect in the random singlet phase: The RSRG 
calculation predicts that all components of the spin and 
string order correlations should decay with the same ex- 
ponents although the underlying XXZ Hamiltonian is not 
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FIG. 12: Exact diagonalization results for the x and z com- 
ponents of the string order at disorder strengths W — 0.5, 
W = 1 and S = 10. 



rotationally invariant. This follows from the fact that the 
ground state of two S-l/2 spins coupled together by an 
interaction of the form 



(36) 



is a rotationally invariant singlet, independently of the 
anisotropy A. So if all the spins really are bound pair- 
wise in singlets, then the decay of different components 
of the correlation functions should be identical. However, 
at finite disorder strength, although the components are 
found to decay with the same exponents, the prefactors 
are different >^> This is due to the fact that for finite disor- 
der strength the strong bonds in the system are not nec- 
essarily surrounded by much weaker bonds, which leads 
to fluctuations in the singlet couplings. As the disorder 
strength is increased these fluctuations should diminish 
and true rotational invariance should be observed. Since 
the string order converges fairly quickly to the expected 
random singlet exponents it is a suitable quantity to use 
to illustrate this crossover behavior. In Fig. Ijl2(l the x 
and z components of the string order are shown for the 
XX chain with disorder parameters W = 0.5, W — 1 and 
S = 10. For flat disorder of strength W = 0.5 the string 
correlation functions already decay with the expected RS 
exponent, but the two components are quite far from each 
other. Increasing the disorder strength to W — 1 the two 
components approach each other, and for a power law 
distribution with <5 = 10 they are within about 10% of 
each other. Increasing S further brings them closer still, 
but it is necessary to use very high numerical precision 
to get reliable data. 

In order to check the decay of the string order away 
from the XX point we again use the SSE method. Here 
we will use chains of length L = 256 and go to sufficiently 
low temperatures to observe T — > converged string cor- 
relations. In Fig. the temperature effects are illus- 
trated for an XX system at disorder strength W = 0.5. 
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FIG. 13: SSE results for the string correlations of an L = 256 
XX system at W = 0.5 calculated at inverse temperatures j3 = 
2" with n = 0, . . . , 15. The (3 = 2 15 results are shown with 
solid circles; the string correlations decrease with increasing 
temperature (decreasing n). . 



In this case, it is possible to obtain T — > converged 
results for all distances. For W close to 1, this would 
require prohibitively low temperatures, but it is still pos- 
sible to obtain well converged results up to distances suf- 
ficiently long for observing the asymptotic RS behavior. 
Fig.^]also illustrates that the string correlations, unlike 
spin-spin correlations, are not symmetric with respect to 
r = L/2 in these periodic systems. From the definition, 
Eq. ©, it is clear that S(r) cannot be symmetric unless 
the total magnetization ^ Sf — 0. This is the case in 
the ground state, where indeed the symmetry is observed. 

In Fig. El low-temperature results are shown at dif- 
ferent W. Here deviations from the RS behavior due to 
temperature effects can be seen for r > 20 when W = 1, 
whereas deviations due to effects of the periodic bound- 
aries (flattening out close to r = L/2) can be seen at 
W = 0.5. In Fig. El we show similar results for the 
XXZ chain for two different combinations of the Ising 
anisotropy A and the disorder strength W. In both cases 
RS behavior can be observed over a significant distance 
range, before temperature or boundary effects become 
visible for r > 50. The very good agreement with the 
RS exponent provides further evidence that the system 
indeed is in the RS phase for any anisotropy and disorder 
strength. 



V. CONCLUSION 

In this paper we have investigated numerically the 
spin-i antiferromagnetic random-exchange XX and XXZ 
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FIG. 14: SSE results for the string correlations of an L = 256 
XX system at different disorder strengths W , calculated at 
the inverse temperatures indicated in the figure. The straight 
line shows the T = RS power law. 
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FIG. 15: SSE results for the string correlations of L = 256 
XXZ system at two combinations of Ising anisotropy A and 
disorder strength W. The line shows the T = RS behavior. 



chains for varying disorder strength. Using exact diag- 
onalization calculations at the XX point and quantum 
Monte Carlo SSE simulations for A > we studied the 
ground state spin-spin and string correlation functions 
for system sizes up to L = 4096 for A = 0. With the 
SSE calculations for A > we went up to L — 256 and 
down to very low temperature, for instance we reached 
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Pmax = 2 17 at the random XXX point for L = 96 and 
disorder strength W = 0.6). We found clear evidences 
for the asymptotically universal behavior of the corre- 
lation functions as predicted by the RSRG analysis of 
Fisheri The main issue of our work presented here con- 
sists in the detailed analysis of the RG flow from the 
pure instable line of XXZ fixed points toward the attrac- 
tive infinite randomness fixed point. Indeed, as we have 
demonstrated, such a flow is controlled by a disorder de- 
pendent length scale £ which diverges as the randomness 
approaches zero^S. In our large scale numerical calcula- 
tions we showed that the spin-spin correlation function is 
very sensitive to such crossover effects whereas the string 
order converges more rapidly to its asymptotic RS value. 
Nevertheless the string order also displays a crossover 
phenomena, visible not in the decay exponents as in the 
spin-spin case but rather in the prefactors. 

The spin-spin correlation function as well as the stiff- 
ness display a clear crossover from the pure behavior 
to the IRFP behavior as predicted by the RSRG. The 
crossover length scale, extracted from numerical data, is 
shown to diverge as £ ~ D"" 1 '. Our estimates for the 
exponent 7 sa 1.0 agrees very well within the error bars 
with the localization length exponent calculated within 
an analytical bosonization approach 6 . However, as the 
bosonization approach is only valid for a disorder that 
is not too strong, our estimates for the crossover length 
scale £(2?) and for the localization length £*(2?) both de- 
viate (in a perfectly similar way) from the predicted be- 
havior described by Ea. (|14f> when the disorder strength 



increases. For any strength of randomness, we found a 
better parameter to describe crossover as well as local- 
ization effects. Indeed, using the variance of the loga- 
rithm of the random couplings, S given by Ea. ll9|) . our 
estimates for the crossover length scale fits in the whole 
range of disorder strengths considered here very well the 
form £((5) - £*(£) oc <T* with $ = It would be 

interesting to check such a 6 dependence of £ or £* also 
for A ^ or 1. The connection between crossover and 
localization effects has been clearly demonstrated here 
and has motivated further studies of the localization in 
ID 30 . 

Whereas the models we have studied are described by 
the IRFP for any strength of the disorder, several dis- 
ordered magnetic systems require a critical value of ran- 
domness to display universal RSP features. For instance, 
gaped systems like the spin-1 chains or spin-i n-legs lad- 
ders are not unstable with respect to the introduction of 
weak disorder and a precise identification of the critical 
disorder T> c might be easier if one considers the diver- 
gence of £ when the disorder strength approaches the 
critical value T> — > T) c . 
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